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We analyze the velocity distribution function of force-free granular gases in the regime of homo- 
geneous cooling when deviations from the Maxwellian distribution may be accounted only by leading 
term in the Sonine polynomial expansion. These are quantified by the magnitude of the coefficient 
a2 of the second term of the expansion. In our study we go beyond the linear approximation for a2 
and observe that there are three different values (three roots) for this coefficient which correspond 
to a scaling solution to the Boltzmann equation. The stability analysis performed showed, however, 
that among these three roots only one corresponds to a stable scaling solution. This is very close to 
02, obtained in previous studies in a linear with respect to 02 approximation. 
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The granular gases, i.e., rarefied systems composed 
of inelastically colliding particles have been of particu- 
lar interest during the last decade (e.g. [0-0 ). Com- 
pared to gases of elastically colliding particles, the dissi- 
pation of energy at inelastic collisions leads to some novel 
phenomena in these systems. One can mention clus- 
tering (e.g. ||]), formation of vortex patterns (e.g. [^), 
etc. Before clustering starts, the granular gas being ini- 
tially homogeneous, keeps for some time its homogeneity, 
although its temperature permanently decreases. This 
regime is called the homogeneous cooling regime (HC). 

In the present study we address the properties of the 
velocity distribution of granular particles in the regime of 
HC, such as the deviation from the Maxwellian distribu- 
tion and the stability of the distribution function. We as- 
sume that the restitution coefhcient e does not depend on 
the impact velocity, i.e. that e = const. The properties of 
the velocity distribution for the system with the impact- 
velocity dependent coefficient of restitution (e.g. g) will 
be addressed elsewhere M. 

It is well known that granular gases in the HC regime 
do not reveal MaxwelHan distribution (^.g. |||,^,^,|j). The 
high- velocity tail is overpopulated [^,||, while the main 
part of the distribution is described by the sum of the 
Maxwellian one and the correction to it, written in terms 
of the Sonine polynomial expansion (e.g. |||,^,|^). Usually 
only the leading, second term, in this expansion is taken 
into account [y|Q,|3l, moreover in previous studies [g[Q] 
only linear analysis with respect to the coefficient a2, 
which refers to this second term has been performed. 
Finding that a2, obtained within the linear approxima- 
tion, is small, the authors of Ref. y,^ conclude a poste- 
riori that the the linear approximation is valid. 

In our approach we also assume that one can restrict 
oneself to the leading term in the Sonine polynomial ex- 
pansion and ignore the other. However we go beyond 
the linear approximation with respect to the coefficient 
02 and perform complete analysis within this level of the 



system description. We observed that there are three dif- 
ferent values of 02 exist which correspond to the scaling 
solution of the Boltzmann equation. The stability analy- 
sis for the velocity distribution function shows, however, 
that only one value of a2 corresponds to a physically ac- 
ceptable stable scaling solution. The stable solution is 
close to the result previously obtained within the linear 
analysis Q]. 

To introduce notations and specify the problem we 
briefly sketch the derivation of the coefficient 02 [3,||, 
in accordance with the approach developed in Ref. [ 3 1 . 

So far we introduce the (time-dependent) temperature 
T(i), and thermal velocity vo{t), which are related to the 
velocity distribution function /(v,i) as 



^nT{t)^Jd^r^fi^r,t)^^nvl{t) 



(1) 



here n is the number density of the granular gas, the par- 
ticles are assumed to be of a unit mass (m = 1), and (|l|) 
is written for 3D-systems. The inelasticity of collisions is 
characterized by the coefficient of the normal restitution 
e (here we consider smooth particles) , which relates after- 
coUisional velocities v*, Vj to the pre-collisional ones, vi, 
V2 as: 



4/2 = V1/2T 2(l+e)(vi2 



e e 



(2) 



where V12 = vi — V2 is the relative velocity, the unit 
vector e = ri2/|ri2| gives the direction of the vector 
ri2 = ri — r2 at the instant of the collision. The time- 
evolution of the velocity distribution function is sub- 
jected to the Enskog-Boltzmann equation, which for the 
force- free case reads |||,^: 

^/(v, t) = 52(o-)o-^ / dv2 / dee(-vi2 • e)|vi2 • e| 



X ^/(vr,t)/(v2**,t)-/(vi,i)/(v2,t) 



(3) 



where a is the diameter of particles, (72(0') — (2— 77)/2(l — 
r;)3 (?7 = i Trncr"^ is packing fraction) denotes the contact 
value of the two-particle correlation function ijl^l , which 
accounts for the increasing collision frequency due to the 
excluded volume effects; 9(x) is the Heaviside function. 
The velocities v^* and Vj* refer to the precoUisional ve- 
locities of the so-called inverse collision, which results 
with vi and V2 as the after-collisional velocities. The 
factor 1/e^ in the gain term appears respectively from 
the Jacobian of the transformation dv^dvj* -^ dvidv2 
and from the relation between the lengths of the coUi- 
sional cyhnders e|v5'2 • e\dt = |vi2 • e\dt ||,p. 

Assuming that the velocity distribution function is of 
a scaling form: 



/(v,i) 



vm 



/(c) 



(4) 



one can show, that the scaling function satisfies the time- 
independent equation 0] : 



fu 



'^i[)f<''^''{f^f) 



with the dimensionless collision integral: 



(5) 



I [f, f)= I dc2 / dee(-ci2 • e)|ci2 • e| 

X {6-v~(cr)/(cr)-.f(ci)/(c2)} (6) 

and with the moments of the dimensionless collision in- 
tegral §: 



A*p = - / dci(?j{j,pj 



while the time-evolution of temperature reads: 
dT/dt^ -(2/3)BT^2 



(7) 



(8) 



where B = B{t) = v^{t)g2{<j)a'^n. 

To proceed we use the Sonine polynomial expansion 
for the velocity distribution function 0,0] 






(9) 



where 0(c) = tt '^/'^ exp(— c^) is the Maxwellian distribu- 
tion and the first few Sonine polynomials read: So{x) = 

1, Si{x) = -x^ + |, S2{x) = 4 - f^ + f ' *5t'=- ^i^^ti- 
plying both sides of Eq. dq) with c^ and integrating over 
dci, we obtain 0]: 



—p {c^} = flp 



(10) 



where integration by parts has been performed and where 
we define 



cPf{c,t)dc. 



(11) 



The odd moments (c^""*"^) are zero, while the even ones, 
/c^"), may be expressed in terms of ak with < fc < n. 
Calculations show that (c^) = |, implying ai = 0, ac- 
cording to the definition of the temperature (nl) (e.g. p|), 
and that (c^) = ^(l + aa). 

Now we assume, that the dissipation is not large, so 
that the deviation from the Maxwellian distribution may 
be accurately described only by the second term in the 
expansion (m with all high-order terms with p > 2 dis- 
carded. Then (nG) is an equation for the coefficient 02. 
Using the above results for (c^) and (c'*) it is easy to 
show that Eq. dl^) converts for p = 2 into identity, while 
for p = 4 it reads: 



5^2 (1 -f 02) - ^4 = 



(12) 



The coefficients fip may be expressed in terms of 02 due 
to the definition (R) and the assumption / = 4>{c)[l + 
025*2(2^)]. Using the properties of the collision integral 
one obtains for fip B: 

^p = -- / dci / dc2 / dee(-ci2 • e)|ci2 • e|0(ci)(/)(c2) 

{l + a2[S2{cl) + S2{cl)]+alS2{cl)S2{cl)}A{cJi + cP2) 

where Aip{ci) = [^(c*) — i^{ci)] denotes change of some 
function ip{ci) in a direct collision. Calculations, similar 
to that, described in 0, yield the following result (some 
detail are given in 0): 



^J'2 



\/2^(l - e') 1 
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0-2 
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1024 



and 



At4 = 4\/2^ {Ti + a2T2 + aJTs} 



with 



Ti 



;(l-^') 



(13) 



(14) 



(15) 



T2 = (1 

128^ 



1 



e^)(69 + 10e^) + -(l + e) 



T3 = — (l + e)H 

^ 64:^ ^ 8192 



(l-e2)(9-30e2) 



The coefficients fi2 and fi4 were provided in Ref. B up to 
terms of the order of 0(02). One obtains the coefficient 
02 in the Sonine polynomial expansion in this approxi- 
mation by substituting (|lj,|lj) into (|lj) and discarding 
in Eqs. (p^|p^) all terms of the order of 0{al): 



a^ 



16(l-e)(l-2e2) 
81-17e + 30e2(l-e) 



(16) 



Calculations including the next order terms 0(0^) in 
the coefficients fj,2 and /i4 show that Eq. (Oh is a cu- 
bic equation, which for physical values of e, < e < 1, 



has three different real roots, as it shown on Fig. ||. 
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FIG. 1. The left hand side of Eq. ^ over a-z for e = 0.8. 
Obviously Eq. [13 has three real solutions. 

Although the cubic equation may be generally solved, 
the resultant expressions for the roots are too cumber- 
some to be written explicitly. However, one of the roots 
(the middle one) is rather small and close to that given by 
Eq. (|l^) , obtained within the linear approximation. This 
suggests the perturbative solution of the cubic equation 
near this root: 



NE 

02 = a^ 



^ _ 1005(1 -e^)-4096r3 ^NE 
6080(1- £2) -4096T2"^ 



(17) 



where we do not write explicitly terms of the order 
O {[02^]'^) and high-order terms. In Fig. g the depen- 
dence of a^^ and of the corresponding improved value 
02 are shown as a function of the restitution coefhcient 
e. As one can see from Fig. the maximal deviation be- 
tween these is less than 10% at small e and decreases as 
e tends to 1. 
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FIG. 2. The second Sonine coefficient 02 as a function of 
the coefficient of restitution e (full line). The dashed line 
shows 02*^ in the first order approximation by van Noiie and 
Ernst [3] according to Eq. (|lq). The approximation ( |l7[ ) is 
shown by circles. 



The other two roots, shown on Fig. || are of the or- 
der of 1 or 10, i.e. are not small. Physically, this means 
that one can not cut the Sonine polynomial expansion 
in this case at the second term and next order terms 
are not negligible to be discarded. Taking into account 
the next order terms, i.e., releasing the assumption that 
Op ~ for p > 2, breaks down the above analysis, since 
the coefficients ^2, M4 occur to be dependent not only 
on a2, but on 03, 04, . . . as well. Thus the occurrence 
of several roots for the 02, found within the above ap- 
proach, which satisfy the conditions required by the scal- 
ing ansatz (|j) does not imply the existence of several 
different scaling solutions. Nevertheless such possibility 
may not be completely excluded. If one assumes that few 
scaling distributions of the velocity may realize, depend- 
ing on the initial conditions at which the HC state has 
been prepared, a natural question arises: Whether the 
particular scaling solution is stable with respect to small 
perturbations, and what is the domain of attraction of 
this particular scaling solution in some parametric space. 
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FIG. 3. The other two solutions for second Sonine coeffi- 
cient 02 of Eq. O over the coefficient of restitution e. 

Certainly, the stability problem is very complicated to 
be solved in general. Therefore, we restrict ourselves to 
the stability analysis of the scaling distribution (Q) where 
the scaling function /(c) has nonzero value of the coef- 
ficient a2, while the other coefficients a^ with p > 2 are 
negligibly small. (For this scaling solution our above re- 
sults for the coefficients /i2, /i4 are valid). Moreover, 
we assume, that small perturbations of the (vanishingly 
small) coefficients a^ with p > 2 do not influence the sta- 
bility of the distribution, and analyze the stability only 
with respect to variation of the coefficient 02 . 

To analyze the stability of the velocity distribution we 
write it in a more general form: 



/(v,t) 



^l(i) 



f{c,t) 



(18) 



which leads, as it easy to show, to the following general- 
ization of Eq. (I) §]:' 



3 



3 + ci 



dci 



/(c, t) + B-'^J{c, t) = i (/, /) (19) 



with the coUisional integral and coefficients fip being now 
time-dependent. The quantities (c^) also depend now on 
time. The temperature, however, evolves still according 

to (D. 

Using / = 0(c) [1 + a2{t)S2{c'^)] and performing es- 
sentially the same manipulations which led before to 
Eq. (|lj), we arrive at the following equation for the co- 
efficient 02 (t): 



d2 - (4/3) B^i2 (1 + 02) + (4/15) 5^4 = 



(20) 



with /X2, Hi still given by ( |l3|Jl4| ), but with the time- 
dependent coefficient 02 (i). Writing the above value B{t) 
as 



B{t) = {87rr'/\c{0)-\{ty/^ 
r,(0)-l=4^l/252(a)a2n^o^/^ 



(21) 
(22) 



where Tc(0) is related to the initial mean-collision time 
at the initial temperature Tq, and u{t) = T{t)/To is the 
reduced temperature, we recast Eq. (po| ) into the form: 



da2 
di 



^^\'/'Fia2) 



15 



(23) 



high accuracy by Eqs. ( p7| , p^ , and with negligibly small 
other coefficients 03, 04, . . . of the Sonine polynomial ex- 
pansion is a stable one with respect to (relatively) small 
perturbations. 

In conclusion, we analyzed the velocity distribution 
function of the granular gas with a constant restitution 
coefficient at the regime of the homogeneous cooling. We 
assume that the deviations from the Maxwellian distri- 
bution may be described using only the leading term in 
the Sonine polynomial expansion, with all other high- 
order terms discarded. In this approach the deviations 
from the Maxwellian distribution are completely char- 
acterized by the magnitude of the coefficient 02 of the 
leading term. We go beyond previous linear theories and 
perform a complete analysis (on the level of the descrip- 
tion chosen), without discarding any nonlinear with re- 
spect to 02 terms. 

Performing the stability analysis of the scaling solution 
of the Enskog-Boltzmann equation we observe, that only 
one value of 02, obtained within our nonlinear analysis 
corresponds to a stable scaling solution. We also report a 
corrections for this value of 02 with respect to the previ- 
ous result of the linear theory. This corrections are small 
(less than 10%) for all values of the restitution coefficient 
e and vanishes as e tends to unity in the elastic limit. 



where t is the reduced time, measured in units of rc(0), 
and where we define a function: 



F{a2) = 5^2(1 + 02) - fJ-i ■ 



(24) 



The form of the function F{a2) for some particular value 
of e is shown on Fig. |l|. This form of F{a2) persists for all 
physical values of the restitution coefficient, < e < 1. 
This has three different roots, ^(03) — 0, i ~ 1,2,3, 
which makes da2/dt vanish yielding the scaling form for 
the solution of the Enskog-Boltzmann equation. The sta- 

. (i) 

bility of the scaling solution, corresponding to Oj re- 
quires for the derivative dF/da2, taken at Oj to be neg- 
ative, since only in this case a small deviation 02 — 02 
from Oj , corresponding to a scaling solution will decay 
with time. As one can see from Fig. |l| only the mid- 
dle root, which corresponds to small values of 02, and is 
close to a^^, predicted by linear theory Q, has negative 
dF/da2, and thus is stable. We also observed that for any 
< e < 1 the point 02 = belongs to the attractive inter- 
val of this stable root. Naturally, this means that initial 
Maxwellian distribution will relax to the non-Maxwellian 
with 02 ~ 0,2^- 

Note that relaxation of any (small) perturbation to this 
value of 02 occurs, as it follows from Eq. (p3|), on the col- 
lision time-scale, i.e., practically "immediately" on the 
time-scale which describes the evolution of the tempera- 
ture. Therefore we conclude, that the scaling solution of 
the Enskog-Boltzmann equation with 02 corresponding 
to the middle root of the function F{a2), given with a 
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